
capture program drop causal_est
program define causal_est
 
 
set seed 1234
global bin=10
global binname=$bin
global bin_new=$bin * 2 
global binname_new=$bin * 2 

global cutoff=300
global data_in=20
global data_out=1180
global range             "RC>=float(${data_in}) & RC<=float(${data_out})"
global range_rv            "incbin${binname}>= float(${data_in}+${bin}) & incbin${binname}<=float(${data_out}+${bin})"
global numerate=20

global zoom_in_c300     =80      //160     
global zoom_out_c300    =470
global range_c300      "incbin${binname}>=float(${zoom_in_c300}) & incbin${binname}<=float(${zoom_out_c300})"


 
matrix ITT_Table =J(2,4,.)  


****************************************************
//(1) keep data
use "${processed_data}\trading_rights_asif.dta" , clear 
rename $myvar y
if $winsor ==1 {
 _pctile  y,  p(1 99)
drop if y<=r(r1) | y>=r(r2) 
}
drop if y==.
keep      if  ${firm}==1 & ${condition}  &  ${year_condition}  & ${range} 
gen tinctobin= RC/${bin}  //registered capital; unit 10K RMB
gen incbinaux=floor(tinctobin)
gen incbin${binname}=incbinaux*${bin}
drop  tinctobin   incbinaux
keep if ${range_c300}

gen RC_2=RC*RC
gen RC_3=RC_2*RC
gen D=0
replace D=1 if incbin${binname}>=float(${cutoff})
gen round50=0
replace round50=1  if mod(incbin${binname},50)==0
gen round100=0
replace round100=1  if  mod(incbin${binname},100)==0  
gen round50X  =round50*RC
gen round50X2 =round50*RC_2
gen round100X =round100*RC
gen round100X2=round100*RC_2
global rounding_q  "round50 round100 round50X  round100X  round50X2 round100X2"


keep y id year incbin* RC*  D  $rounding_q  cic_adj2 
save "${causal_data_dir}/data_for_estimation_${myvar}.dta",replace



****************************************************
//Step 2: Main Regression with TE estimates  

use "${causal_data_dir}/data_for_estimation_${myvar}.dta", clear 
qui tab cic_adj2, gen (sector)	
reg y  D   RC   RC_2   $rounding_q  sector*    if (incbin${binname}<float(${graph_low}) ) | (incbin${binname}>float(${graph_high}) ), cluster(incbin${binname})			  
predict y_ct  
//for bootstrap starts
predict res, res
replace res=0 if incbin${binname}>=float(${graph_low}) & incbin${binname}<=float(${graph_high})  
gen y_fit=y_ct if incbin${binname}<float(${graph_low})  | incbin${binname}>float(${graph_high})
replace y_fit=y if incbin${binname}>=float(${graph_low}) & incbin${binname}<=float(${graph_high})
save "${causal_data_dir}/data_for_bootstrap_${myvar}.dta",replace	
//for bootstrap	ends   
collapse (mean) y  y_ct  D=D , by(incbin${binname}) 
merge 1:1 incbin${binname} using "${density_data_dir}\density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta"
drop  numbin${binname}_noround numbin${binname}_counter_noround round50 Z100  _merge
replace numbin${binname}_counter=round(numbin${binname}_counter)
sort incbin${binname} 
qui keep if incbin${binname}>=float(${graph_low}) & incbin${binname}<=float(${graph_high})

qui su y_ct [fweight=numbin${binname}_counter]
matrix ITT_Table[1,2]=round(r(mean),   0.0001)

 
gen numbin${bin}_complier=.
replace  numbin${bin}_complier=numbin${bin}- numbin${bin}_counter if incbin${bin}>=float(${cutoff})  
replace  numbin${bin}_complier=numbin${bin}_counter-numbin${bin}  if incbin${bin}<float(${cutoff})
replace  numbin${bin}_complier=0 if numbin${bin}_complier<float(0) 

gen y_complier=.
replace y_complier= (numbin${binname}*y - numbin${binname}_counter*y_ct) / (numbin${binname}_complier)      if incbin${binname}>=float(${cutoff})
replace y_complier= (numbin${binname}_counter*y_ct - numbin${binname}*y) / (numbin${binname}_complier)      if incbin${binname}<float(${cutoff})
if $trim ==1 {
replace y_complier=0 if y_complier<0
replace y_complier=1 if y_complier>1 & y_complier~=. 
}

replace  numbin${bin}_complier=round(numbin${bin}_complier)


su numbin${binname}_complier                         if incbin${binname}<float(${cutoff})
global sum_cp_ct=r(sum)  
gen test_cp_ct=y_complier* numbin${binname}_complier if incbin${binname}<float(${cutoff})
su test_cp_ct
global E_cp_ct=r(sum)/${sum_cp_ct}
drop test_cp_ct


su numbin${binname}_complier                      if incbin${binname}>=float(${cutoff})
global sum_cp=r(sum)  
gen test_cp=y_complier* numbin${binname}_complier if incbin${binname}>=float(${cutoff})
su test_cp 
global E_cp_ob=r(sum)/${sum_cp}
drop test_cp

global ITT_cp=${E_cp_ob} - ${E_cp_ct}
matrix ITT_Table[1, 3]=${ITT_cp}

gen y_cp_mean=.
replace y_cp_mean= $E_cp_ct if incbin${binname}<float(${cutoff})
replace y_cp_mean =$E_cp_ob if incbin${binname}>=float(${cutoff})

if $draw_fig ==1{
drop y y_ct D numbin${binname} numbin${binname}_counter 
global bin_new=$bin * 2  
gen tinctobin= incbin${binname}/${bin_new} //draw figure at bin=20
gen incbinaux=floor(tinctobin)
gen incbin${bin_new}=incbinaux*${bin_new}
drop  tinctobin   incbinaux

bysort incbin${bin_new}: egen  numbin${bin_new}_complier= sum(numbin${binname}_complier)
gen temp1=numbin${binname}_complie*y_complier 
bysort incbin${bin_new}: egen temp2=sum(temp1)
gen temp3=temp2/numbin${bin_new}_complier
rename temp3   y_complier_new
drop temp*
bysort incbin${bin_new}: gen dup=_n
drop if dup>1
drop dup y_complier incbin${binname} numbin${binname}_complier
				 
if $high <=10 &  $bin_new ==20 {
set obs 6
replace y_cp_mean=y_cp_mean[5] if _n==6  
replace incbin${bin_new}=320   if _n==6
}
#delimit;
twoway
(scatter  y_complier_new  incbin${bin_new}  [fweight= numbin${bin_new}_complier],   mcolor(black   black)  msymbol(o)   msize(large large))
(line  y_cp_mean  incbin${bin_new} if incbin${bin_new}<float(${cutoff}),  sort clcolor(black)   lwidth(thick) lpattern(dash))
(line  y_cp_mean  incbin${bin_new} if incbin${bin_new}>=float(${cutoff}), sort clcolor(black)   lwidth(thick) lpattern(solid))
,
legend( order (2 "Treated Average of Compliers"  3 "Counterfactual Average of Compliers" )    size(vsmall))
xline(${cutoff}, lcolor(black) lwidth(thin))
ytitle(${myvar_name}, size(small)  margin(medium))  
xtitle(Registered Capital (unit:10,000 RMB) in ${year_name}, size(small) margin(medium) ) xtitle(, alignment(top)) 									
ylabel(#5) 
graphregion(fcolor(white)  ifcolor(white) color(white) icolor(white))
title("${myvar_name}", justification(center) color(black) size(medsmall)) 
subtitle( "for Compliers of ${firm_name} firms in ${location} ${year_note}"  "${cutoff_note}",  justification(center) color(black) size(medidum))
note("Note: bin=${bin_new}; data: ASIF.",size(small) )
;
graph export   "${causal_graph_dir}/fitfig_causal_compliers_${cutoff}_${firm_name}_${name}_${year_name}_${myvar}.pdf", replace;
#delimit cr 	
}		


******************************************************************
/*step 2: bootstrap */

if $run_bootstrap ==1 {

forvalues ii=1 (1) $nboot {
display `ii'

use "${causal_data_dir}/data_for_bootstrap_${myvar}.dta", clear
keep y_fit res id year incbin* RC*  D  $rounding_q cic_adj2 
qui su year, detail
global rows=r(N)
gen draw_j=ceil(uniform()*${rows})
gen resample=res[draw_j]  
gen y=y_fit+resample
drop y_fit  res  draw_j  resample   


qui tab cic_adj2, gen (sector)	
reg y  D   RC   RC_2   $rounding_q  sector*    if (incbin${binname}<float(${graph_low}) ) | (incbin${binname}>float(${graph_high})) , cluster(incbin${binname})
predict y_ct 
collapse (mean) y=y   y_ct=y_ct   D=D , by(incbin${binname}) 
merge 1:1 incbin${binname} using "${density_data_dir}\density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta"
drop  numbin${binname}_noround numbin${binname}_counter_noround round50 Z100  _merge
replace numbin${binname}_counter=round(numbin${binname}_counter)
sort incbin${binname} 
qui keep if incbin${binname}>=float(${graph_low}) & incbin${binname}<=float(${graph_high})

gen numbin${bin}_complier=.
replace  numbin${bin}_complier=numbin${bin}- numbin${bin}_counter if incbin${bin}>=float(${cutoff})
replace  numbin${bin}_complier=numbin${bin}_counter-numbin${bin}  if incbin${bin}<float(${cutoff})
replace  numbin${bin}_complier=0 if numbin${bin}_complier<float(0)  

gen y_complier=.
replace y_complier= (numbin${binname}*y - numbin${binname}_counter*y_ct) / (numbin${binname}_complier)      if incbin${binname}>=float(${cutoff})
replace y_complier= (numbin${binname}_counter*y_ct - numbin${binname}*y) / (numbin${binname}_complier)      if incbin${binname}<float(${cutoff})
if $trim ==1{
replace y_complier=0 if y_complier<0
replace y_complier=1 if y_complier>1 & y_complier~=. 
}

replace  numbin${bin}_complier=round( numbin${bin}_complier)

su numbin${binname}_complier                         if incbin${binname}<float(${cutoff})
global sum_cp_ct=r(sum)  
gen test_cp_ct=y_complier* numbin${binname}_complier if incbin${binname}<float(${cutoff})
su test_cp_ct
gen E_cp_ct_`ii'=r(sum)/${sum_cp_ct}
drop test_cp_ct

su numbin${binname}_complier                      if incbin${binname}>=float(${cutoff})
global sum_cp=r(sum)  
gen test_cp=y_complier* numbin${binname}_complier if incbin${binname}>=float(${cutoff})
su test_cp 
gen E_cp_ob_`ii'=r(sum)/${sum_cp}
drop test_cp

gen ITT_cp_`ii'=E_cp_ob_`ii' - E_cp_ct_`ii'
drop  E_cp_ob_`ii'
keep if _n==1 
keep   ITT_cp_`ii' 
gen year=${year_name}
qui save  "${causal_temp_dir}/data_bt_est_${myvar}_boot`ii'.dta",replace
}


/* sort bootstrap results*/
use "${causal_temp_dir}/data_bt_est_${myvar}_boot1.dta", clear 
foreach ii of numlist 2 (1) $nboot {
merge 1:1 year using "${causal_temp_dir}/data_bt_est_${myvar}_boot`ii'.dta"
drop _merge 
}
save  "${causal_bt_dir}/BT_est_${density}_reg${reg_degree}_${fe_status}_${cutoff}_${firm_name}_${name}_${year_name}_${myvar}_degree${degree}_${low}_${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta",replace


preserve 
keep ITT_cp_*
xpose, clear 
qui su v1 , detail
matrix ITT_Table[2, 3]=round(r(sd), 0.0001)
matrix ITT_Table[2, 4]= round(ITT_Table[1,3] / ITT_Table[2,3],0.0001)
restore


foreach ii of numlist 1(1) $nboot {
qui rm "${causal_temp_dir}/data_bt_est_${myvar}_boot`ii'.dta"
}


svmat  ITT_Table
keep   ITT_Table* 
rename ITT_Table1  varname 
rename ITT_Table2  DepMean
rename ITT_Table3  TE
rename ITT_Table4  t 
tostring varname, replace
replace varname="${myvar}" if _n==1
replace varname=""         if _n==2
keep if _n<=2
save  "${causal_ITT_dir}/ITT_${density}_reg${reg_degree}_${fe_status}_${cutoff}_${firm_name}_${name}_${year_name}_${myvar}_degree${degree}_${low}_${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta", replace
qui rm "${causal_data_dir}/data_for_bootstrap_${myvar}.dta"
}

qui rm "${causal_data_dir}/data_for_estimation_${myvar}.dta"

end




